library(lmtest)
library(sandwich)
library(doBy)


load(file="data/sf_2.Rdata")##full data set 2008-2015, all suspected crimes


#CONSTRUCT LAGGED DV
means<-summaryBy(found.weapon~dayno, data=s, FUN=mean, na.rm=T)##get daily hit rates
means$lag.hit.new<-NA
means<-means[order(means$dayno, decreasing=F ),]
head(means)

for(i in 2:length(means$dayno)){
	
	means$lag.hit.new[i]<-means$found.weapon.mean[i-1]
	
}


head(means)


print(dim(s))
s<-merge(s, means, by="dayno", all.x=T)
print(dim(s))

s<-s[order(s$dayno),]

print(identical(s$dayno, s$dayno[order(s$dayno)]))
stopifnot(identical(s$dayno, s$dayno[order(s$dayno)]))



##TABLE 1
## Estimate discontinuity at moment of intervention using all weapon stops, 2008-2015

se1<-rep(NA, 8)
names.models<-c("diff.means","diff.means+controls","linear","linear+controls","squared","squared+controls","cubic","cubic+controls")
names(se1)<-names.models

##diff in means
m1<-lm(found.weapon~period, data=s)
se1[1]<-summary(m1)$coefficients["periodpost-memo","Std. Error"]
m1a<-lm(found.weapon~period+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=s)
se1[2]<-summary(m1a)$coefficients["periodpost-memo","Std. Error"]

##local linear
m2<-lm(found.weapon~period+distance+period*distance, data=s)
se1[3]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]
m2a<-lm(found.weapon~period+distance+period*distance+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=s)
se1[4]<-summary(m2a)$coefficients["periodpost-memo","Std. Error"]

##second order poly
m3<-lm(found.weapon~period+distance+period*distance+I(distance^2)+period*I(distance^2), data=s)
se1[5]<-summary(m3)$coefficients["periodpost-memo","Std. Error"]
m3a<-lm(found.weapon~period+distance+period*distance+I(distance^2)+period*I(distance^2)+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=s)
se1[6]<-summary(m3a)$coefficients["periodpost-memo","Std. Error"]


##cubic
m4<-lm(found.weapon~period+distance+period*distance+I(distance^2)+period*I(distance^2)+I(distance^3)+period*I(distance^3), data=s)
se1[7]<-summary(m4)$coefficients["periodpost-memo","Std. Error"]

m4a<-lm(found.weapon~period+distance+period*distance+I(distance^2)+period*I(distance^2)+I(distance^3)+period*I(distance^3)+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=s)
se1[8]<-summary(m4a)$coefficients["periodpost-memo","Std. Error"]


models2<-list(m1, m1a, m2, m2a, m3, m3a, m4, m4a)

print("models done")

save(se1,file= "output/se_std_alldata.Rdata")##save conventional SEs
save(models2,file= "output/std_models_alldata.Rdata")##save conventional SEs

vc<-list(NA)

##estimate HAC SEs

print("computing HAC SEs")


for(i in 1:length(models2)){#reset to 1 if replicating
vci<-vcovHAC(models2[[i]])
vc[[i]]<-vci
save(vc, file="output/all_data_HAC_vc.Rdata")
print(i)
print(vc[[i]][2,2])
}

